Over the past few decades, the lattice Boltzmann (LB) method has been proven to be an efficient alternative approach in many computational fluid dynamics (CFD) areas, e.g., hydrodynamic systems,^[1,2] multiphase and multicomponent fluids,^[3–7] and porous media flow.^[8,9] With the rapid development of multi-core supercomputer and parallel computation, the LB method is gaining more attention from CFD researchers due to its inherent parallelism.^[10–12] Historically, the LB method arose from the pioneering idea of lattice-gas automata,^[13] however more recently we tend to link it with the kinetic equation, the Bhatnagar–Gross–Krook (BGK) Boltzmann equation,^[14] which leads to the reconstruction of LB theory.^[15–21]

Although the newly constructed LB theory has abandoned the original lattice-gas theory, interest in reproducing the classical lattice BGK (LBGK) equation inherited from lattice-gas automata, which has been proven numerically stable, is still considerable. To achieve this, numerous approaches have been proposed, e.g., Strang splitting,^[15] Maxwell iteration,^[18] Taylor expansion,^[20,22] etc. However, the approximation and truncation employed in these approaches make the reproduction unconvincing. Taking the most popular Taylor expansion scheme as an example, the LB equation is dealt with as a characteristic integral of the BGK Boltzmann equation with constant collision term. To fix the truncation error introduced by the constant collision term assumption, the scheme expands the left-hand-side integral of the BGK Boltzmann equation in a Taylor series while recovering the Navier–Stokes equations with a Chapmann–Enskog expansion.^[22] Employing a 2nd-degree Taylor expansion, the LBGK equation can be properly recovered. But the truncated Taylor series introduce uncertainty. And a further research report^[16] asserts that by employing the Taylor expansion scheme to recover the Navier–Stokes energy equation, the viscosity in viscous heat dissipation term is inconsistent with the momentum equation. Hence He et al.^[16] employed a trapezoidal rule to improve the approximate accuracy of the collision term to avoid Taylor expansion. Recently, Boesch and Karlin^[19] used an Euler–Maclaurin integral to analytically integrate the collision term, eliminating the truncation error of the trapezoidal rule. Compared with the trapezoidal rule, the result of the Euler–Maclaurin integral (exact LB equation as called in the paper) is theoretically more rigorous. The problem is that the elegant but extremely complicated Euler–Maclaurin integral cloaks the important effect of equilibrium distribution. It makes the derivation in Ref. [19] be founded on a quite rough approximation of equilibrium distribution, which leads Boesch and Karlin to the conclusion that the classical LBGK equation cannot be recovered from the analytical characteristic integral of the BGK Boltzman equation. In particular, when the viscosity is minimal, the relaxation time in the LBGK equation would be greater than 1 tending to 2, i.e., over-relaxation, and meanwhile it is always under 1 in the exact LB equation. Thus Boesch and Karlin^[19] asserted that the over-relaxation in the LBGK equation is beyond the kinetic theory of the continuous-time BGK Boltzmann equation. Then an essential question arises: can the LBGK equation be exactly reproduced based on the BGK Boltzmann equation?

In this paper, we propose an LB theory to link the LB method and the BGK-Boltzmann equation, designated as analytical characteristics integral (ACI) LB theory. ACI theory follows the philosophy of Refs. [16] and [19] to construct the LB equation, i.e., improving the accuracy of the collision term integral to eliminate the truncation error and inconsistency in the Taylor expansion scheme. Through analyzing the BGK Boltzmann equation, ACI theory asserts that the BGK-Boltzmann equation is accurate enough to describe the particle kinetics behind the Navier–Stokes continuum and momentum equations, and the LB equation can be constructed via analytically integrating the BGK-Boltzmann equation along characteristics with approximated temporal evolution of equilibrium distribution. ACI theory employs differential equation solving skills to achieve the analytical integral. Compared with the Euler–Maclaurin integral employed in Ref. [19], differential equation solving skill is mathematically more concise with a clear physical figure. Meanwhile the rigorous characteristic integral avoids the truncation error of trapezoidal rule in Ref. [16]. In ACI theory, the temporal evolution of equilibrium distribution along characteristics is the kernel of constructing the LB equation, determining the computational stability and accuracy of LB equations. It extends the generality of BGK-Boltzmann kinetic theory discussed in Ref. [19]. To demonstrate it, we recover several popular LB equations and numerically analyze them. Our derivation also indicates that the characteristic parameter of an LB equation is collision number instead of the traditionally assumed relaxation time, which is merely a reflection of the temporal evolution of equilibrium distribution along characteristics, and the over-relaxation time problem is a manifestation of the evolution model behind the LBGK equation under high collision number, supported by the kinetic theory of the continuous-time BGK-Boltzmann equation.

This paper is organized as follows. In Section 2, we propose the ACI LB theory and recover popular LB equations. In Section 3, we demonstrate the discretization of the LB equation under ACI LB theory. In Section 4, we numerically analyze the derived LB equations with several simple benchmarks. Section 5 concludes the paper.

The BGK-Boltzmann equation is a result of linear approximation on the collision term of the Boltzmann equation^[14]

The BGK-Boltzmann equation is a kinetic equation describing the transient interaction between particles. In LB theory, we use it to depict the particle kinetics behind the Navier–Stokes continuum and momentum equations. It can be mathematically validated with classical Chapman–Enskog (CE) expansion.^[23] CE expansion is a multiple scale technique, proposed to solve the Boltzmann equation at first.^[23] The key concept of CE expansion is formulating the flow phenomena with different time scales. It decomposes the partial {derivatives} with respect to t and r into different scales

Now, we substitute the expansions in Eq. (4) and Eq. (5) into the BGK-Boltzmann equation in Eq. (1) and collect the terms with the same order of ε to generate new equations on different scales. The first two scaling equations of ε read

It should be noted that this recovery of the Navier–Stokes equations do not involve Taylor expansion, avoiding the truncation error. The derivation shows that the BGK-Boltzmann equation is accurate enough to describe the particle kinetics behind the Navier–Stokes continuum and momentum equations. Assuming the BGK-Boltzmann equation is the exact kinetic equation for these two equations, all we need to do is construct the LB equation based on the BGK-Boltzmann equation. Here we follow the philosophy of the Boesh–Karlin approach,^[19] i.e., constructing the LB equation through analytically integrating the BGK-Boltzmann equation along characteristics. Actually, with simple equation transformation,^[20] all terms in the BGK-Boltzmann equation can be directly integrated except the equilibrium term (Maxwell–Boltzmann distribution). We start with transforming the left side of the BGK-Boltzmann equation, treating it as the temporal derivative along characteristic line ξ,

The integration in Eq. (15) has not been finished, where the integral of g(t′) still remains unsolved. Since the calculation of g(t′) involves macroscopic density ρ and velocity u (the temperature T would be dealt with as a field constant in practice), which evolve nonlinearly and need to be calculated by the LB equation, it is impossible to get its analytical distribution along the characteristics. The only solution is using approximation. Since the BGK-Boltzmann equation is assumed to be exact for the particle kinetics behind the Navier–Stokes continuum and momentum equations, and the LB equation is its analytical integration along characteristics, then the accuracy of the LB equation solely depends on the design of g(t′) evolution along characteristics, i.e., the aforementioned temporal evolution of equilibrium distribution. In other words, the g(t′) model forges the final form of the LB equation. It greatly extends the kinetic theory of the BGK-Boltzmann equation discussed in the Boesh–Karlin approach,^[19] which we will demonstrate in detail in the following. In case the evolution of the BGK-Boltzmann equation departs from the physical process too much, we restrict the value of g(t′) at t′ = 0 to gⁿ, ensuring that the transient equation of particle interaction with respect to (r, ξ, t) is an exact BGK-Boltzmann equation. We call this LB theory the analytical characteristic integral (ACI) LB theory.

We start with typical treatment for an unknown distribution, i.e., assuming g(t′) is constant over [0,Δt], namely, the steady assumption (SA) model,

Now the question is what is the g(t′) model behind the LBGK equation. As we assume that the LB equation is an analytical characteristic integral of the BGK-Boltzmann equation and the g(t′) model ensures that the transient equation of particle interaction with respect to (r, ξ, t) is an exact BGK-Boltzmann equation, then the g(t′) model is not difficult to deduce, namely, the evolutionary correction of deviation (ECD) model,

Summarizing all LB equations derived in this section, the SA, Boesh–Karlin, ECD, DCD and He–Chen–Dong LB equations, they all are analytical characteristic integrals of the BGK-Boltzmann equation. The difference is in their designs of g(t′) along the characteristics. In practice, Boesh–Karlin and He–Chen–Dong LB equations equal SA and ECD LB equations, respectively. All these LB equations are implemented with a unified form

We close this section with a few remarks.

(I) In this section, we propose the ACI LB theory. In ACI LB theory, we take the BGK-Boltzmann equation as the exact kinetic equation behind the Navier–Stokes continuum and momentum equations and construct the LB equation through rigorously integrating the BGK-Boltzmann equation along characteristics. The BGK-Boltzmann equation can be analytically integrated except a Maxwell–Boltzmann distribution, which evolves nonlinearly and needs to be calculated by the LB equation. Thus the design of g(t′), the temporal evolution of equilibrium distribution along characteristics, decides the accuracy of the LB equation. In case the whole evolution of the BGK-Boltzmann equation departs from the physical process too much, we restrict the value of g(t′) at t′ = 0 to gⁿ, ensuring that the transient equation of particle interaction with respect to (r, ξ, t) is an exact BGK-Boltzmann equation.

(II) In ACI LB theory, there are two important parameters, the collision time λ, depicting the time required for the distribution to reach equilibrium through collision, and the step time Δt, denoting the time interval of the LB equation. Their ratio, the collision number Δt/λ, describes the times of the distribution reaching equilibrium state. It reflects the intensity of particle interaction in the time span of the LB equation, independent from detailed g(t′) models. In contrast, the traditionally assumed characteristic parameter of the LB equation, relaxation time, is merely a reflection of the g(t′) model. Thus the actual characteristic parameter of the LB equation should be collision number Δt/λ, which is supported by Eq. (32).

(III) ACI LB theory is a physically general theory. By tuning the g(t′) model, ACI LB theory can easily reproduce the SA, Boesh–Karlin, ECD, DCD, and He–Chen–Dong LB equations which are developed under different theories. It allows the comparison of LB equations to avoid ambiguous theoretical analysis of their derivations, leads to a simple investigation of their designs of g(t′) model. Another direct implication of this generality is that the form of LB equation is by no means limited to Eq. (31) since the flexibility of g(t′) model is far greater than the ones presented in this section.

(IV) ACI LB theory is a mathematically rigorous but concise theory. It avoids the viscosity inconsistency^[16] introduced by truncating Taylor series in Taylor expansion schemes.^[20,22] As an analytical characteristic integral approach like the Boesch–Karlin approach, ACI theory eliminates the uncertainty of the trapezoidal rule in the He–Chen–Dong approach^[16] (see Ref. [19] for more detailed discussions). Compared with the Boesch–Karlin approach,^[19] the integrating skill of ACI theory is mathematically more concise with a clear physical figure. And the asserted concept that the g(t′) model forges the LB equation significantly extends the understanding of BGK-Boltzmann kinetic theory in the Boesch–Karlin approach. ACI theory also avoids the uncertainty of collision and streaming splitting operation in the Strang splitting approach.^[15]

(V) ACI LB theory is sensitive. The differences among LB equations will clearly manifest in their corresponding g(t′) models. For example, the g(t′) models corresponding to the Boesh–Karlin and He–Chen–Dong LB equations are different from their implemented equations, i.e., SA and ECD LB equations.

In the former section, we have demonstrated the derivation of the LB equation from the BGK-Boltzmann equation. But the whole argument is based on a velocity-continuous equation, which should be discretized before implementation. In ACI LB theory, the discretization of an LB equation in velocity space is independent of the construction of the LB equation; it is directly performed on the BGK-Boltzmann equation. The discretization includes discretizing the velocity space of particle distribution and constructing the equilibrium distribution based on the discrete-velocity model. Technically speaking, the discrete equilibrium distribution model determines the final recovery form of hydrodynamic equations, such as compressible or incompressible equations, conduct equation, and so on, but this is out of this paper’s discussion. Here we offer a detailed demonstration and illustrate the framework of ACI theory, i.e., the discretization of velocity space on the BGK-Boltzmann equation, determining the practical LB algorithm, is independent of the detailed LB equation. We take a classical discretization as the detailed demonstration, small Mach-number approximation (SMA) approach, proposed by Ref. [20]. For the sake of simplicity but without losing generality, we employ a 2-dimensional Maxwell–Boltzmann distribution to construct the classical D2Q9 model. The derivation can be easily extended to 3 dimensions and new models constructed.

The key idea of velocity discretization in ACI LB theory is to ensure that the assumed macroscopic hydrodynamic equations can still be recovered by the discretized BGK-Boltzmann equation. To achieve this, the SMA approach employs keeping the relative velocity moment integrals of discretized equilibrium distribution consistent with the Maxwell–Boltzmann distribution. For instance, the following velocity moment integrals in Eq. (33) have been used in Section 2 to recover the Navier–Stokes continuum and momentum equations:

To construct the discrete equilibrium distribution, the SMA approach firstly decomposes the exponent part in the Maxwell–Boltzmann distribution: one for the weight function, the other for the distribution with respect to macro-velocity,

For the convenience of calculating weights, the truncated small velocity expansion (or small-Mach-number approximation) is employed,^[20] where the terms above 2nd macro-velocity order are neglected,

After small-Mach-number approximation, the target continuous equilibrium distribution is g′ instead of the original Maxwell–Boltzmann distribution g. Right now, g′ in Eq. (35) bears a strong resemblance to the classical D2Q9 equilibrium distribution, and all that remains to be accomplished is calculating the weight of discrete velocity. As we discussed before, the discrete equilibrium model f^eq should keep the velocity moment integrals required in Eq. (33) consistent with g′,

The derived D2Q9 model is directly constructed on the BGK-Boltzmann equation, independent of the detailed LB equation. Thus it is applicable to all g(t′) models.

Within ACI LB theory, all popular LB equations are analytical characteristic integrals of the BGK-Boltzmann equation, identified by their designs of g(t′) model, so then all we need to analyze is their designs of g(t′) model. ACI LB theory also proposes a new characteristic parameter of a LB equation, collision number Δt/λ, depicting the particle-interacting intensity in the time span of the LB equation. Unlike the traditionally assumed characteristic parameter, relaxation time 1/τ, which is merely a reflection of g(t′) model, the collision number is independent of detailed g(t′) models. In this section, we analyze the numerical performance of derived LB equations under different collision numbers. In numerical computation, the Boesh–Karlin and He–Chen–Dong LB equations would be replaced by the SA and ECD LB equations, ignoring the implicitness as their proposers suggested.^[16,19] We start with the comparison of relaxation time curves along collision number as all derived LB equations share a unified equation form in Eq. (31), only distinguished in relaxation time formulas. Then two benchmarks, unsteady Couette flow and lid-driven cavity flow, will be calculated by SA and ECD LB equations, employing the derived D2Q9 model. The results of SA and ECD will be compared to analyze their accuracy under different collision numbers and validate the effect of g(t′) model. In the calculation of benchmarks, the DCD approximation will be abandoned due to its instability.

4.1. Relaxation time

As Eq. (32) shows, for the same calculation, i.e., the same collision number Δt/λ, the proposed g(t′) models, SA, DCD and ECD, only differ in relaxation time 1/τ. Figure 1 illustrates the relaxation time profiles of SA, DCD and ECD models with respect to collision number. For Δt/λ → 0, SA, DCD and ECD approach each other. As the collision number increases, all three relaxation times increase with different rates, which leads to 1/τ_DCD > 1/τ_ECD > 1/τ_SA. When the collision number goes to infinity, the SA and ECD relaxation times converge to 1 and 2, respectively, while that of DCD increases unlimitedly, which causes the instability of the DCD model. Figure 1 indicates that SA, DCD and ECD LB equations would have the same numerical performance under low collision number but are distinguished from each other under high collision number. This conclusion is supported by the design of g(t′). Once collision number is small enough, which indicates small collision change, the particle distribution f(t′) along characteristics would not depart from fⁿ too much. Then the modification parts in ECD model Eq. (23) and DCD model Eq. (25) would both approach 0,

Hence the g(t′) models of ECD and DCD can be approximated as which is exactly the SA model. When the collision change increases as the collision number increases, the modification parts in Eq. (54) are not neglectable. Thus SA, ECD, and DCD models would show different results. It is worth noting that under large collision change, the SA and DCD models are theoretically unsuitable as SA takes the equilibrium distribution constant and DCD takes the collision change rate constant, which contradict the physical process. What is more, figure 1 also indicates that the over-relaxation time problem is merely a manifestation of the g(t′) model under high collision number, irrelevant to particle kinetics.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Relaxation time vs. collision number. The solid, dashed, and dashdot lines are the ECD, DCD, and SA results, respectively. Note that the result of DCD is truncated by the limitation of relaxation time.

4.2. Unsteady Couette flow

Unsteady Couette flow, also known as transient plane Couette flow, is a typical numerical benchmark. It evolves in a straight channel, which infinitely extends in the x direction. The walls of the channel are parallel to the x axis and defined by the equations y = 0 for the lower wall and y = L for the upper wall. The flow is stationary at the beginning, and driven by the upper wall with a constant velocity after that. Its equation can be described as

with the following initial and boundary conditions in the range 0 ≤ y ≤ L: where ν is kinematic viscosity, U is the upper wall’s driven velocity, and L is the width of the channel. Unsteady Couette flow can be solved analytically,

The periodic boundary is implemented in the x direction to simulate infinite extension. For the upper and lower walls, the regularized boundary^[24] is applied. Actually, in unsteady Couette flow, the Zou–He boundary^[25] and the Inamuro boundary^[26] would share the same macro-result with regularized boundary, though the particle distributions f may differ. Five cases are designed to discuss the SA and ECD models (see Table 1). Each case has two states: low collision number (LCN) and high collision number (HCN). Figure 2 shows the SA and ECD velocity profiles of unsteady Couette flow at t = 1.00 × 10⁻³ s. Figure 3 illustrates the error evolution in time. The error formula is defined as

where u is the computing velocity profile, and u_ana is the analytical velocity profile. For the sake of space saving, only the error evolutions of case 1, case 3 and case 5 are plotted.

Fig. 2.

Figure Option ViewDownloadNew Window

Fig. 2. The velocity profiles of unsteady Couette flow at t = 1.00 × 10−3 s: (a) low collision number; (b) high collision number. The symbols ■, ▲, ⋆, • and ♦ are the results of case 1, case 2, case 3, case 4 and case 5, respectively. Under low collision number in (a), the models agree well with each other, and share a similar profile with the analytical solution. Under high collision number in (b), ECD stays consistent with the analytical solution while SA fails. The detailed descriptions of cases are listed in Table 1.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. Velocity error evolution of unsteady Couette flow: (a) case 1, (b) case 3, (c) case 5. ECD-HCN, ECD-LCN, SA-HCN, SA-LCN represent the results of the ECD model with low collision number, ECD model with high collision number, SA model with low collision number, and SA model with high collision number, respectively.

Table 1.

Table 1.

Table 1. Configurations of unsteady Couette flow demonstration cases. For all cases, the channel width equals L = 0.1; m, and the driven velocity of the upper wall is U = 1.0; m/s. . Case CN Grid Time step Kinematic Collision type number viscosity number case 1 LCN 201 2.50×10−4 1.00×10−3 3.33×10−1 HCN 2.50×10−6 3.33×101 case 2 LCN 101 1.00×10−4 1.00×10−2 3.33×10−1 HCN 1.00×10−6 3.33×101 case 3 LCN 101 2.00×10−5 5.00×10−2 3.33×10−1 HCN 2.00×10−7 3.33×101 case 4 LCN 101 1.00×10−5 1.00×10−1 3.33×10−1 HCN 1.00×10−7 3.33×101 case 5 LCN 101 1.00×10−6 1.00×100 3.33×10−1 HCN 1.00×10−8 3.33×101

Table 1. Configurations of unsteady Couette flow demonstration cases. For all cases, the channel width equals L = 0.1; m, and the driven velocity of the upper wall is U = 1.0; m/s. .

The velocity profiles of the SA LB equation in Fig. 2 show that at LCN, SA results are close to the analytical solutions, while at HCN, they evidently deviate, manifesting magnified viscosities. In contrast, ECD results always stay consistent with analytical solutions regardless of collision number. The deviation of the SA LB equation under HCN is explainable, due to the intensive collision change; the equilibrium distribution would be greatly changed, which leads to the failure of the constant equilibrium distribution assumption in the SA model. And the results of the ECD LB equation indicate that the modification part in the ECD model significantly improves its numerical performance under HCN. Comparing the results of SA and ECD LB equations, it is clearly shown that g(t′) model determines the accuracy of the LB equation. The consistent numerical performance of SA and ECD LB equations under different collision numbers also validates the fact that the collision number is a good characteristic parameter to depict LB equations. It reflects the collision intensity in the time span of the LB equation, independent of detailed g(t′) model.

4.3. Lid-driven cavity flow

Lid-driven cavity flow is a classical constant-property benchmark for fluid computation, because of its simple geometry and complicated flow behaviors, especially the corner flow phenomena. Lid-driven cavity flow is a 2-D case characterized by Reynolds number (Re), which develops in a square cavity with side length L. The fluid would keep stationary at the beginning and be driven by the top lid with constant velocity after that. Except the top lid, all other walls remain at rest. This case has no analytic solution. As a convention, the results in Ref. [27] are taken as a calibration on velocity profiles. As the definition of collision number shows, with small kinematic viscosity, the collision number is inherently large. It is very difficult to simulate a high Re-number cavity flow under low collision number without mesh explosion. Fortunately, our aim is to evaluate g(t′) models under collision number instead of Re number. A small Re-number case with different collision numbers will fulfill the task. Hence the case with smallest Re number in Ref. [27], Re = 100, is calculated for analysis.

The simulation is configured as shown in Table 2. The cavity sides and corners employ the regularized boundary.^[24] Figure 4 illustrates the velocity profiles along the vertical and horizontal lines through the cavity center. The models keep consistent with their performance in Subsection 4.2.: under LCN, SA and ECD they share a close profile; under HCN, ECD sticks with the LCN profile while SA fails. To give a full picture of the calculations, the vorticity contours are plotted in Fig. 5, which can be quantified by post-processing software compared to the streamlined case. The values of vorticity along contours are 0, ± 0.5, ± 1.0, ± 2.0, ± 3.0, 4.0, ±5.0. As Fig. 5 illustrates, under LCN, SA (see Fig. 5(b)) and ECD (see Fig. 5(a)) have the same vorticity contour; under HCN, although the contour of ECD in Fig. 5(c) is a bit dispersed, the shape agrees well with the LCN result while for SA (see Fig. 5(d)) it is quite different.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (a) Velocity profiles of Ux along the vertical line through the cavity center. (b) Velocity profiles of Uy along the horizontal line through the cavity center. The line labeled Ghia et al. is the calibration from Ref. [27].

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Vorticity contours of the cavity flow. The values of vorticity along contours are 0, ± 0.5, ± 1.0, ± 2.0, ± 3.0, ± 4.0, 5.0.

Table 2.

Table 2.

Table 2. Configurations of lid-driven cavity flow. The cavity side length equals 1.0 m, and the driven-lid velocity is 1.0 m/s. . CN Grid Time step Kinematic Collision number viscosity number LCN 257 1.00×10−3 1.00×10−2 5.09 × 10−1 HCN 1.00×10−6 5.09×10+2

Table 2. Configurations of lid-driven cavity flow. The cavity side length equals 1.0 m, and the driven-lid velocity is 1.0 m/s. .

The results of the SA model illuminate the fact that, under LCN, as the collision change is minimal during a timestep, the constant equilibrium distribution assumption in the SA LB equation can properly approximate its temporal evolution, while under HCN, due to the great collision change, it seriously deviates from the physical temporal evolution. The accurate result of the ECD model under HCN indicates that the modification part in the ECD model could accurately approximate the change of equilibrium distribution caused by intensive collision. The different numerical performances of SA and ECD LB equations under HCN confirm the effect of g(t′) model on the accuracy of the LB equation. Compared with unsteady Couette flow in Subsection 4.2., the governing equations of lid-driven cavity flow, i.e., full incompressible-flow Navier–Stokes equations, are more complicated and are close to a real CFD problem. Thus the consistence of models’ numerical performances between unsteady Couette and lid-driven cavity flow validates the universality of our assertion on the accuracy of SA and ECD LB equations under different collision numbers. Hence, we can conclude that, for a fluid calculation, SA and ECD are both available under LCN while ECD is a better choice under HCN.

In this paper, we propose a mathematically rigorous but physically general LB theory, ACI LB theory. In ACI LB theory, we assume that the BGK-Boltzmann equation is an exact kinetic equation behind the Navier–Stokes continuum and momentum equations and construct the LB equation by analytically integrating the BGK-Boltzmann equation along characteristics. Since the integral of g(t′) along characteristics cannot be analytically solved, the approximation is employed. To ensure that the evolution of the LB equation does not depart from the physical process too much, we restrict the value of g(t′) at t′ = 0 to gⁿ so that the transient equation of particle interaction with respect to (r, ξ, t) is an exact BGK-Boltzmann equation. In ACI LB theory, the g(t′) model is the determinant of LB equation accuracy, and collision number is the characteristic parameter of the equation, depicting the particle-interaction intensity in the time span. As a demonstration, we recover the approximating models behind popular LB equations, e.g., SA, Boesh–Karlin, ECD, DCD and He–Chen–Dong LB equations, and numerically analyze them, where Boesh–Karlin and He–Chen–Dong equations are replaced by SA and ECD since they are equivalent, respectively, after removing the implicitness of gⁿ⁺¹.

To conclude, we highlight the following conclusions drawn from our work: (i) An LB equation can be constructed through analytically integrating the BGK-Boltzmann equation along characteristics with approximated temporal evolution of equilibrium distribution g(t′), avoiding the approximation and truncation in classical schemes.^{[15,16,20,22]} (ii) The kinetic theory depicted by the BGK-Boltzmann equation is quite general, far beyond the discussion in the Boesh–Karlin scheme.^[19] By tuning the g(t′) model, most popular LB equations can be recovered as analytical characteristic integrals of the BGK-Boltzmann equation, including the questioned LBGK equation. (iii) An LB equation is identified by its g(t′) model. The approximation accuracy of the g(t′) model determines the LB equation accuracy. Among the g(t′) models presented, ECD is a better choice for numerical simulations. (iv) The characteristic parameter of an LB equation is collision number Δt/λ, instead of the traditional relaxation time 1/τ, which is merely a reflection of g(t′) model. It depicts the particle-interaction intensity in the time span of the LB equation, independent of detailed g(t′) model. And the over-relaxation time problem is merely a manifestation of ECD model under high collision number.